/*  slmkhmlt.cpp  by K.Tsuru  */
// function ID = 220    DRADIX, BRADIX
/*********************************************************************
SLong and SInteger classes
 z = x*y
It provides the Karatsuba's multiplication using recursive calls.
Knuth
x*y=(R^n*x1+x0)*(R^n*y1+y0)
  = {R^(2n)+R^n}*(x1*y1) + R^n*(x1-x0)*(y0-y1) +(R^n +1)*(x0*y0)
It is desirable that x and y have the identical size.
**********************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

void SLong::KHHMult(const SLong& x, const SLong& y, SLong& z){
    if( (&z == &x) || (&z == &y)){
        z.SetError(z.SYNTAX_ERR, "KHHMult", 220);
    }
    if( !x.Sign(220) || !y.Sign(220) ){
        z.SetZero();    return;
    }
    uint n = x.FFTMaxArraySize()/2;
    uint xfig = ceilpow2(x.Head()+1), yfig = ceilpow2(y.Head()+1), k;
    if( (xfig <= n) && (yfig <= n) ){
        z = x*y;    // FFT
        return;
    }

    SLong x0(x.Type(), 0), x1(x.Type(), 0);
    SLong y0(x.Type(), 0), y1(x.Type(), 0), u(x.Type(), 0);

// If xfig != yfig, divide x or y into two parts.
    if(xfig > yfig){
    // (R^n*x1+x0)*y = R^n*(x1*y) + x0*y, where n = xfig/2
        SLDivHalf(x, x0, x1);
        KHHMult(x0, y, z);
        KHHMult(x1, y, u);
        z += u.MultPowRdx(xfig/2);
        return;
    } else if(xfig < yfig){
        SLDivHalf(y, y0, y1);
        KHHMult(y0, x, z);
        KHHMult(y1, x, u);
        z += u.MultPowRdx(yfig/2);
        return;
    }
// Here xfig == yfig is satisfied.
    SLDivHalf(x, x0, x1);
    SLDivHalf(y, y0, y1);
    k = yfig/2;
    KHHMult(x0, y0, u);
    z = u;
    z += u.MultPowRdx(k);
    KHHMult( x1 - x0, y0 - y1, u);
    z += u.MultPowRdx(k);
    KHHMult(x1, y1, u);
    z += u.MultPowRdx(k);
    z += u.MultPowRdx(k);
    z.KaratsubaHHMultUsedTimes++; // counter
    return;
}